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ABSTRACT 

A doubly curved element for a shell of revolution which has arbitrar- 
ily curved meridians is developed and analyzed. Menicaaaen curvature is 
calculated using a highly accurate polynomial approximation. The dis- 
placement functions selected satisfy interelement compatibility and 
contain all the lower modes of a complete set of straining eu cee : 
Non-straining modes corresponding to rigid body motions are introduced 
into the final stiffness matrix for any conceivable rigid body motions. 

The direct stiffness method was used to construct a stress and 
strain analysis program and the results of the analysis of a few problems 
are compared to classical solutions to establish the integrity of the 


element. 
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I, INTRODUCTION 


A. GENERAL DESCRIPTION 

The concept of dividing complex structures into an assemblage of 
individual structural components or elements is axiomatic to the struc- 
tural engineer. However, the ability to execute a breaking down and 
reassembling process has been dependent upon the availability of an- 
alytical tools. Not until the present decade were developments in high- 
speed digital computers sufficient to make possible the fundamental 
approach to problems of structural analysis known as the finite element 
method. 

The finite element method is essentially the analysis of an idealized 
structural system composed of a finite number of elements interconnected 
ata finite number of joints or nodal points. It is the finite characteristic 
of structural connectivity which distinguishes the finite element problem 


from problems of continuum mechanics. 


B. HISTORICAL BACKGROUND 

The finite element method of analysis is applicable to one, two and 
three dimensional structures. Background discussion will be limited to 
the two dimensional case, although recently there has been considerable 
research activity on three dimensional structures. 

Various shapes of finite elements have been employed in plane stress 
and plane strain analysis. Early research was focused on the development 


of flat rectangular and triangular finite elements. An abundance of papers, 
v 


technical reports, and theses were produced; a complete bibliography 
would be beyond the scope of this paper. It is sufficient to say however, 
that requirements for efficient flat plate finite elements have been formu- 
lated and tested. Considerable success was reported in Refs. 1 and 2 
using these elements to approximate shells of free form. 

The development of curved finite elements, a considerably more 
difficult problem, has progressed at a much slower pace. Some early 
success was achieved in the development of truncated conical ring 
elements. However, problem formulation using these elements has not 
been generalized. Truncated conical elements have produced excellent 
results for particular problems but have also resulted in unsatisfactory 
analysis of similar problems. 

A ring element with arbitrarily curved meridians developed by R. E. 
Jones and D. R. Strome [Ref. 3] produced acceptable results for the case 
of shells of revolution under pressure loads. 

In 1968, G. Cantin [Ref. 4] successfully developed a curved finite 
element for cylindrical shells. This curved element was developed using 
6 nodal coordinates per joint (u,, Vi, Wye Wee « Wins Wye n) and implicitly 
contained a complete description of rigid body motions. This element was 
tested and found to produce excellent results for a broad class of cylin- 
drical shell and flat plate problems for both membrane and bending 


behavior. (The flat element is obtained by letting R~®&.) 


C. WEJECTIVE 

The objective of the author's research is to develop the stiffness 
matrix of a four nodal point finite element for a shell of revolution which 
has arbitrarily curved meridians, hereafter to be called a curvilinear shell 


finite element. 


TL LOE LLL Se Eevi Nave Tr @ab 
A. GENERAL 
Finite element analysis of an elastic continuum logically separates 
into four phases: (1) structural idealization, (2) evaluation of the element 
stiffness properties, (3) the evaluation of the force-displacement be- 
havior of the overall element assemblage, and (4) the evaluation of 


element stress and strain. 


B. STRUCTURAL IDEALIZATION 

Structural idealization is the division of the actual elastic continuum 
into an assemblage of geometrically compatible structural elements which 
are interconnected at a diScrete number of nodal points, through which 
the forces are transmitted. 

The discretization of the continuum is a physical process, subject 
only to conditions of geometric compatibility, rather than an approxima- 
tion of a mathematical nature. The idealized structure, then, is con- 
structed of elements which possess the material properties of the actual 
COMME. 

A basic reduction is necessary to make the idealized structure 
mathematically tractable. This reduction, based on an a priori selection 
of possible deformations of the elastic continuum, constrains the 
elements to deform in only a certain number of predictable shapes. These 
constraints called displacement functions uniquely define displacements 
and the state of strain within the element in terms of nodal displacements. 


These functions are the keystone of a successful idealization. 
10 


C,. EVALUATION OF ELEMENT STIFFNESS PROPERTIES 

The evaluation of the stiffness properties of the element is the crit- 
ical phase of the finite element method. As the evaluation of the stiffness 
properties of a curvilinear shell finite element is the object of the 
author's research, procedures will be discussed in following sections. 
For the present, it is sufficient to say that the element stiffness matrix 
relates the nodal forces and nodal displacements. This relation takes 
the form 

B= Lk lot (oj!) 

where 

{F} = nodal force vector 

{§} = nodal point displacement vector 

[k] = element stiffness matrix 
D, EVALUATION OF THE FORCE-DISPLACEMENT BEHAVIOR OF THE 

IDEALIZED STRUCTURE 

The force-displacement behavior of the overall element assemblage 
is obtained by the systematic superposition of individual element stiff- 
ness matrices and load vectors. This technique, called the "direct 


stiffness method," is well documented in the literature. 


E, EVALUATION OF ELEMENT STRESS AND STRAIN 
When nodal displacements are known,element strain is directly 
obtained using classical shell theory strain-displacement relationships. 


Element stress is then calculated using constitutive laws. 


Hil 


III, CONSTRUCTION OF THE STIFFNESS MATRIX 


A, GENERAL 

The stiffness matrix for a particular finite elment is the relationship 
between a Set of generalized ae olacements and the corresponding 
generalized forces, 

Standard matrix notation, square brackets denoting rectangular 
matrices, curly brackets representing column matrices, superscript (-1) 
denoting inversion, and Superscript T denoting transposition, will be 


used, 


B. GEOMETRIC DERIVATION 

The geometric description of the curvilinear shell finite element 
must contain sufficient parameters to guarantee that an assemblage of 
elements will form a smooth surface whose slope is continuous and whose 
dimensions at element boundaries match those of the actual continuum. 
Furthermore, the geometric description must be sufficiently general to 
describe shells.of revolution which have meridional inflection of reverse 
curvature as well as shells of degenerate geometry such as cylindrical 
shells and flat plates. An extension of the descriptive procedure used 
in Ref. 3 was found to satisfy the above requirements. 

Geometrical symbols for the curvilinear shell finite element are 
defined in Figure 1. Nondimensional meridional arc length is denoted 


by 


Ws 





NOTE: ALL DIMENSIONS BASED 
ON SHELL MIDSURFACE 


Fig. | DEFINITION OF GEOMETRICAL SYMBOLS 
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ae: 
5 = L (31) 


and nondimensional parallel arc length is denoted by 


n= - (3.42) 


Through the use of a second order polynomial in € for the cosine of the 
angle ¢, 
Ces te Clg c ce: exe (Sa 


and the values of 


pt) r (2) uaa (1) g@) 2.6 (see Figure 1) 


expressions for completely defining element geometry were directly 


formulated. The radius of the element is given by 


rey = (1) + 2 {cos g¢d& = r(0) +2(c,&E +406" + cab?) 


.o] 


(3 725 
and the principal radii of curvature are 
R, = h a £sin ¢ (3.5) 


d¢/dé —_ co, + 2656 


R= (3 . 6) 


sin @ 
where 
Cc, = cos g (2) (3.7) 
6, = 2 CES al eno g(2) - cos of] (3.8) 
c. = 3 [-2 baal + cos g'1) + cos gf) J]. (3.9) 
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It is worthy of note that all values required for these expressions 
are found from engineering drawings or models. 

For particular cases, such as conical shells and cylindrical shells 
where geometry is defined by analytical expressions, geometric 


approximation is not required, 


C. DISPLACEMENT FUNCTIONS 


The displacement functions 


u(é ,n) 


u} (3. 8G) 


v(é ,n) 


w (é ,n) 
must satisfy three critical requirements. 

1. They must insure that continuity of the structure can be auto- 
matically maintained across common boundaries of adjacent elements. 

2. They must insure that all the lower modes of a complete set of 
straining modes will be present in the stiffness matrix, 

3. They must contain an exact description of all possible rigid 
body motions of the element. 

It can be shown [Ref. 5] that satisfaction of the first two requirements 
will guarantee that the strain energy of the "assemblage of elements" 
approximation will represent, for specified loads,a lower bound to the 
Strain energy of the actual continuum. Satisfaction of the third 
requirement is necessary to expose hidden constraints opposing rigid 


body motion, 


lis 


mee Application*to"tite CurvilinearShell Fim te flemens 


The polynomial expansion: 


weyy) =a, err" + axe "deh nets ¥ (3..11lu) 
vig.) = as&n tag& + ay + ag (3.1 1 
w(§,7) = ag an + a6°n* a aii6°n 1 ai2é° (3. l1w) 


2.3 2 2 

oa ais § ” +a,46° of <i ais nN + aes 
3 2 

ae aiy En +ae5N ay aS as Azac g 
S 2 

42 Gael, +a + Aon Ta 


Of Meresecompact!y 

{u}= LP] {a,} (3a) 
satisfies the requirements of interelement compatibility and completeness 
of strain field. Requirement three, that an exact description of all 
possible rigid body motions of the element be represented, is not 


satisfied, however, as will be shown in the following section. 


D alti LD BODY MODES 

It has been established in Ref. 6, that conditions of equilibrium 
are not satisfied unless the displacement functions contain an exact 
description of rigid body motions of the element. Moreover, it has been 
established in Ref. 7 that inclusion of rigid body motion is essential and 
cannot be represented by independent displacement components. Rigid 
body motion, then, must be introduced without compromise to interelement 


compatibility requirements. 
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A general method for including rigid body motion is developed in 
Ref. 6. 
1. Application to the Curvilinear Shell Finite Element 
Consider a point P, specified by position vector p, on the 
element in the system of reference shown in Figure 2. If the element is 
submitted to general translation then, 


Om 6.1 + 65) + 65k (3.13) 


in the system of reference (x, y, z). The corresponding scalar compo- 


nents of the displacement field are; 


u cos ¢ cos 9 cos ¢ sin @ - sin ¢ 6, 
V =] - sin é@ cos 6 0 5. (3.14) 
w sin @ cos @ sing sin 6 cos ¢ 5s 


or more compactly 


{u} = Lee (5, } 


When submitted to rigid body rotation of small amplitude 

B = Bri + Bol + Bok (3.15) 
the displacement vector field of the element in reference system (x, y, Z) 
is 


6, = BX p (Simliey) 


The corresponding scalar expression is 


u - sin 9 (z cos ¢+rsin ¢) cos @(zcos ¢+rsin g) Of; {p, 

vy =| - Z cos 8 -zsin @ 1 Ma fe 

Ww - sin 9(z sin g@-rcos ¢) cos 6@(z sin ¢g-rcos ¢) Of] [Bs 
(L7) 


le 





ae Se WN 


Ce a 
Fig 2 DEFINITION OF COORDINATE SYSTEMS 


or more compactly 


fu} 7 [Ror { B; } 


Combining expressions (3.14) and (3.17) gives: 


u) ' ory 

52 
M 7 Rorans ' ror 8s (3.18) 

By 

F ] Be 

Bs 

or 

fu) = CR] {a} (3.19) 


Matrix [R]is given explicitly in Table l. 


E. OLASSICAL SHELL THEORY 

In the introduction it was emphasized that the difference between 
finite element problems and problems in continuum mechanics is based 
on a physical approximation only. Classical theories for shells from 
continuum mechanics are, therefore, applicable in the finite element 
approach. 

The consitutive law for an isotropic, homogenous, thin elastic shell 


can be written: 


Ny mw ee) omg Oo 7: 

N ae (CG) € 

n ee n 

N 0 oa Gh my € 

én = = En ‘ (3,20) 
M, 0 o> =D. 0 m 

M 0 (ODED, x 

n n } 

M pap eg eo O° x 

bn | (én, 


where 


:Kog = v K,; Ka = ope 


Et 
Di agi) pe VDL: Da = # O-v) D 


or more compactly 
{N} = [E] {e} (Sez) 

Young's modulus E-and Poisson's ratio vy, according to definition 
are constants. The thickness t, however, is allowed to vary in the 
meridional direction. In this thesis, shells of uniform thickness and 
shells of linearly varying thickness are considered .? 

Numerous theories exist for curvilinear shells. The venerable 
theories have been subject to years of critical inspection; their merits 
and limitations are well known. Several theories investigated would 
have been acceptable. Kraus' formulation [Ref. 8] was selected for the 
Pe liewite reasons: 

1. The theory is consistent. That is, the theory holds for the 
limiting cases of cylindrical shell and flat plate elements. This quality 
is of basic importance since mesh size reduction to the SE where the 
element is nearly flat is necessary in many finite element problems. 

2. The theory is presented in a form that can be directly applied 
to this studye 


3. The theory allows strain free modes for rigid body motions. 


1 Any thickness which varies as a function of § could be considered 
by changing only one card in the computer program. 


vAU 


The strain displacement relations are: 














le l 
=r u 
“¢ £ 9& : Ry 
cos leeds sin 
Sn r ar on 6 
13 oo ee 
SE ar On £ as ri as ! 
j 
=! 1 9 1 oR : ae a 
ee oR, Se ” OR MBE se” 
is 2 
cos @¢ pag nem: ” 
*n r Ry ar on anr* 38 
_ os 6 oO. 
r & o& 
l Q a ee eee 2 ar tO | w 
“tn| | ark, 8n ox 86 or d& ras a& a7 
2. 3 
rakegayn 
6S .22) 
or more compactly 
i (Diem) (3.23) 


All classical theories investigated by the author contained inherent 
Singularities, Kraus' formulation was no exception. It can be seen that 
several terms in equation (3.22) become infinite whenever ris equal to 
zero. This theoretical inadequacy has been avoided by specifying either 
a small hole or a small rigid insert with appropriate boundary conditions, 
at the pole point. The problem solved in Chapter 5, Section B, 2 


demonstrates this method. 


Z| 


F, NODAL COORDINATES 

For computational convenience, physical nodal coordinates are 
transformed to a generalized coordinate system. This transformation 
facilitates the expression of displacements, elastic characteristics, and 
strain energy in a system free from geometric complications. The 
following set of nodal coordinates was selected: 

oy = CU, My, Wi, Mie Wie Wis ) c i = 1s 

(3.24) 

where Sand @are the physical coordinates of any point on the element. 
Then using (3.1), (3.2), and (3.12) the transformation matrix [T] is 
CONnswuCteuss UCh that 

fume = fry {a,}. (3.25) 

The transformation matrix [T] is ASheinetier and is easily 
inverted, 


Matrix [T] is given explicitly in Table 2. 


G. STIFFNESS MATRIX 

The strain energy in the element is: 

Us = #{f iN} le}rasdtan (3.26) 
Substituting for stress {N} and strain {¢€} one gets: 

Views $ ‘cml Lam Yoel [P4 (E1(PAre z 9497 ) (tT? {ua 7) 


where 


Pai TE (3.28) 


Ze 


The stiffness matrix is then: 
[kK] = (r*)* ([f CP*] (6) (pede ot 868m (171 


or 

ee) Nae) (3.29) 
where [T]and [k] are both (24x24) matrices. 

Stiffness matrix [K]may contain some rigid body modes depending 
on problem type. However, by subjecting [K] to eigenvalue analysis, 
individual rigid body modes, if present, can be identified. 

When [K] does not contain all six components of rigid body motion, 
it must be modified accordingly. The method developed in Ref. 6 is 


applied as shown below: 
' Uy 
fu,} = E ! R | --5- (3.30) 
J u, 


R T 
@here {u } = (6,,8,,6), Br Boe B, ? . (tl) is a (24x24) 


identity matrix, and [R] is the (24x6) matrix developed in Section D. 
Matrix LK] is then expanded to include all six rigid body modes: 

I 
~~~, LK | JI} R | = | ---se---- eee (3.31) 

rh RD | 


where | (RI [kK] [RI] | is a (6x6) matrix, 


F {k] 


| 
Py ee 


pea al a_i (3.32) 
0 RIKI [RIK] [R] ay 


Solving for 


Zo 


fu} = [tr"txitr2 } [txitri] {v,} (3.33) 


However, if the stiffness matrix [K] already contains some of the rigid 


body modes, then matrix | CR} (K]{R] | will be singular, These 


modes have been previously identified and are merely suppressed in 


{ue } and in the corresponding column of [R]. Matrix | (RI we [RI | 


then becomes nonsingular. Solving (3.33) and substituting into the first 
(3.32) leads to the following expression: 
{ F } = C ies [ Ck]tR] |[ CRI" [K](R] y [trv [K] | ) { uy } 
(3.34) 


The stiffness matrix modified to include rigid body modes is then: 
a +1 1k 
[K* =| kK | - | k]tR) ][ Rv" Ck] (R1)~ | fRV Ck) 


Stiffness matrix [K*] is strain free for any rigid body motion. 


H, DEFINITION OF MATRICES 
In the following sections of this thesis, the stiffness matrix [K*] 


derived above is referred to as CSFE. The matrix | tr] [K](R] | 


developed in (3.31) is hereafter referred to as K22. For the general case 
where all six rigid body modes are absent, K22 is a (6x6) matrix. In the 
case of a cylindrical shell segment where four rigid body modes are 
absent, K22 is (4x4), and in the case of a conical segment where five 


rigid modes are absent, K22 is (5x5), 
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I. LOAD VEGEORS 

Loads and moments acting on the element surface must be repre- 
sented by vectors which are consistent with displacements and rotations 
in the generalized coordinate system. 

A consistent load vector then, is defined as a vector made up of 
fictitious quantities which would give, after an inner product with the 
nodal point displacement, the same work as the real load. 

1. Concentrated Load Vectors 


For concentrated loads of Ps : - , Fy acting at a point (s, 6) 


C 


of the element, the work done is 


W =F,u+F v + Fow 
2 n C 


or more compactly 


wel}? {3 


Substituting for {U} from (3.12) and (3.24) gives, 


w= {o)* Le] Lr] {} 


The consistent concentrated load vector is then: 


{erono} = C7" {r.} 6.35 


2. Pressure Load Vectors 


The work done by a pressure load p is: 


w = | oe als (ana) 
A 


then if 


ae 


fat=(0 0 p)! 


direct substitution of known quantities into (3.36) gives: 


w = || fat’ (8) sar dédn = Cow [ffa¥"[ P][rt |r d&dn ) {us} 


The consistent pressure load vector is then: 


=. tal{[P ]° (q} rdédn. (3.37) 


Table 1. Transformation Matrix [R] Between Nodal Coordinates 


Bae 


and Polynomial Coefficients 


Res, 2 


Ru ,2 


Ri 3 
0 


Ris 3 


Ris,3 


Rig ,3 


Ryo ,3 


Ra 3 


a 


Roe ,3 


Ri ,4 
Bene 


Ris 4 


Ris 4 
Rig 4 
Ri7 4 
Ris 4 
Rie 4 


Reo,4 


Res ,4 


Rea 4 


Ris 5 
0 

Riss 
Rig 5 
Riz 5 
Ries 
Riss 
Rae 5 
Ra js 
Rees 
Ress 


Raa 5 


oe ee Sa. 


wv 


8,6 


eS aS 2] ©} 


Ria,6 


oo 2 }. 2 oS 


Roars 


Co . (Ss 





* Commas in this table are used as Separators, no partial differentia- 


tion is implied. 


i 


Rs 2 
‘ 


Re a 
i 





a r (1) 


2 
Rov + Ro 2+Ri 1 Co +> 
f a 


Zo 





Rio > 
R 

11 : 2 
Ry F 4 


Rye 2 


Rigo 
s 


Ris ,1 
Ris 2 


Rig 3 


Rig 4 = 
a 


Ris 5 
a 
Ria 1 
Ris 2 
Ris 6 
Ris 1 
a 


Ris 2 


Rig ,3 = 


(2) 


(cy + 265)0 


= Ro. + Roig (C, + C, + Ca) + z 


= COs g (2) COs @ 

= cos g (2) sin @ 

= ihe: g (2) 

- sin g'?) sin & r(?) 


(2) 


= sin ¢ COS & rq) 
= - sing 

= cos @ 

= ; (+) 

= sin g (2) COS & 


sin g (2) SiO! 


cos gl) 


Zi 


Ris 4 = Cos 6) sina rr) 
é 
Ris s = - cos 63 cos a rl) 
‘ 
Rig,1 = Ro ,1 Ria 2 
Rig 2 = -Ro a Ria i 
Rig ,3 Rio, 3 
Rig 4 = ra Rig 5 
Mie se oe Ria 2 Rio s 
R 
= 14,1 
Rv 1. ae one aN 
r 
Ris 1 
17,2 ~  +(4) 
Riga = Ris 1 
Rigs = Rig 2 
f a 
R 
14,1 
R =—- 
181 Ryo a (2) 
r 
R 
16,1 
r 
Rie,4 = Rig 2 Ria, 
Rie,6 = = Ria Rie 4 
Rio 1 ~ Raa Ria 2 
Migee) = gees aan 
Rig,s — Rios 
Rie 4 a Ria Riis 
Rigs ~ Ria 2 Riis 
Ree 1 Ria 3 
Rao 2 = Rya,2 
Rao 4 i Re 4 Ria 2 
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Rya 1 Rs,s 
Ryas2 Ra,s 
R R 

a1 Ria,2 


- Raji Risa 


a 


ie : a 
Res a 


me Rya)2 Rg ,4 


ot 


Transformation Matrix |T| Between Nodal Coordinates 


Mable 2: 


(Columns 1 - 12) 


and Polynomial Coefficients 


Cul. 


T, 


oi 
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om 
Tg 


emus 
OT, 
6T; 


or. 
am, 
9T 


cum 


2 4) 
Reus 
lumns 1 

(Co 


d) 
(Continue 


Table 2: 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


Ty 


Ts 


ie 


Tz 


oT 


T, 


T4 


T, 


Tg 


ip 


Tp 
Tz 

Tp 

Tp 

oT. 

rage 

T4 


Zee 
OT, 
aie 


i 
on, 
6T, 


Dae 


V4 
oT 


ot 


Sil 
Sr 


Tg 


De 


7; 


Ty 


qT, 


Ty 


27 


oe 


ig 
2T. 


ate 


where: 


il 
T iach! - = _ 
: ari) = art?) 
PS at co a 1. - Sa 
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IV. NUMERICAL CONSIDERATIONS AND PROGRAM INFORMATION 


A. GENERAL 

All numerical computations were accomplished on the IBM 360/67 
digital Pattie! using FORTRAN 4 language. The entire program was 
coded in double precision arithmetic. An overlay scheme under control 


of the linkage editor was used to conserve core storage. 


B, QUADRATURE 

It was evident that the integrands of (3.27) and (3.37) could not be 
manipulated into closed form. Since the above expressions were in 
nondimensional form, lower and upper limits of integration zero and one 
respectively, it was advantageous to use GausSian quadrature. An 
additional advantage of uSing Gaussian forms is that a polynomial of 
degree 2n-1 can be exactly described by n points [9]. For polynomials 
of unknown degree, maximum accuracy is achieved by selecting a large 
number of points. Additional accuracy, however, is obtained at the 
expense of computation time. Investigation by the author showed that 
Gauss Legendre-—8 point quadrature provided accuracy within the range 
of double precision round off error for shells whose geometry could be 
described by analytic expressions. For shells of more complicated 
geometry, required accuracy is obtained when standard finite element 


mesh size reduction is accomplished, 
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C. LIMITING CASES 

In this thesis, cylindrical shells, truncated cones, and flat plates 
have been considered as special cases. When a Shell is specified as 
one of the above forms, geometric approximation is precluded, 

Exact solutions for these shells provided a powerful tool for the 
analysis of the stiffness properties of the general case of the curvilinear 


shell finite element. 


D,. CASE WHERE ¢ IS EQUAL TO ZERO 

Singularities exist in the geometric approximation whenever the 
element is constructed in such a way that a point on its Surface is 
locates at ¢ equal to zero (see Figure 1). Discussion of the effects of 


this singularity is contained in Chapter 5, Section B, 2. 


E. INITIAL CONSTRUCTION OF THE STIFFNESS MATRIX 

The author originally constructed the stiffness matrix by constructing 
all matrices in (3.27) and then integrating the entire expression. 
Accurate results were achieved. However, the construction of the stiff- 
ness matrix for a single element required 500,000 bytes of computer stor- 
age and approximately 20 minutes of computation time. This was 
obviously unsatisfactory for other than pedagogical application. 

Fortunately, identical results were obtained by individual stiffness 
matrix element computation and by taking advantage of symmetry. This 
was done in the remarkably short time of 15 seconds uSing only 143,000 


bytes of computer storage. 


S116 


CSFE was constructed and evaluated in a separate program prior to 


being incorporated in a direct stiffness program. 


F, DIRECTION STIFFNESS PROGRAM 

Numerical solution of engineering problems was accomplished 
through the use of a stress analysis program which was based on the 
direct stiffness method. The total stiffness matrix and the load vectors 
were assembled uSing programming techniques furnished by Professor G. 
Cantin. This program gave numerical values for displacements, moments, 
stresses, and strains at each nodal point of the mesh. A flow diagram of 
the computer program is presented in Figure 3. It can be seen that the 
program is functionally divided into two parts. In the CSFE portion, the 
artntes smatrix and its associated pressure load, strain-displacement, 
and stress-strain relationships are computed for individual elements. 
Results of these computations are stored on disk units, In the stress 
analysis portion, individual element information is retrieved from disk 
and assembled in a manner consistent with the direct stiffness method. 

The author's primary purpose for construction of the stress analysis 
portion of the program was to verify the stiffness properties of the 
curvilinear shell finite element. It was therefore, designed to handle 
problems of limited size. Specifically, the program is limited to a 
maximum of 36 elements, 49 nodal points, or 294 equations. Coding 
methods used, however, are general to allow for future program size 


expansion. The computer program is found in Appendix C. 
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V. DISCUSSION 


A, INTEGRITY OF THE CURVILINEAR SHELL PINITE ELEMENT 

Confidence in the integrity of CSFE was established by comparing 
stiffness properties of limiting cases (flat plate and cylindrical shell) 
with properties found in the literature. Element by element comparison 
with values given by Bogner, et. al. [Ref. 10] revealed exact agreement 
for the case of the flate plate. In the case of the cylindrical shell 
element similar comparison was made with values given by Cantin [Ref. 
1]; once again exact agreement existed. 

Confidence in the integrity of the stress analysis program was 
established by direct comparison of problem solutions using CSFE with 
solutions found in the literature. Results of comparison are discussed 


in the following sections. 


B, NUMERICAL SOLUTION OF CLASSICAL PROBLEMS 

The classical problems analyzed in this section were selected to 
demonstrate the integrity of the element, and to exhibit both the 
flexibility and limitations associated with CSFE. 

1, Pinehed Cylinder 

The pinched cylinder illustrated in Figure 4 has been analyzed 

by Bogner, etal (Ref. 11] using a (48x48) element and by Cantin 
[Ref. 6] using a (24x24) element. Table 3 is a comparison of results 
obtained using CSFE with those found in the literature. The obvious 


agreement of results confirms the integrity of CSFE. Moreover, itis 
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evident that CSFE has the desirable characteristic of rapid convergence 


to the true solution (-0.1139) given by Cantin. 


No. of Bogner et al Cantin CSFE 
E@s.. (48 x 48) (28 x 28) (24 x 24) 
Mesh Displace. Mesh Displace. Mesh _  Displace. 
(in.) 

54 2CZ -0.0931 Ze -0.0932 
108 2X2 -0.0808 

150 4x4 ~0.1126 4x4 - ,1128 
180 2x4 -0.1098 

294 6x6 -0.1137 6x6 -0.1138 
486 8x8 ~-0.1139 

TOMS 10xl0 -0.1139 





Table 3. Displacement under load P for pinched cylinder 


2. spherical Cap 


The spherical cap illustrated in Figure 5 demonstrates the 
general case of CSFE. This problem which has an exact solution given 
by Timoshenko [Ref. 12] is useful for demonstrating the capabilities and 
limitations of the new element. Figure 6 contains a plot of meridional 


bending moments .* 


“ The signs of Timoshenko's results have been changed to agree with 
the author's convention. 
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It is worthy of note that the pole of the cap was located at a singular 
point (¢ = 0). In order to fulfill connectivity requirements, a small rigid 
polar "cap" was assumed to exist. Boundary conditions were applied to 
allow this "cap" freedom of displacement in the z-direction at zero slope. 
Zero slope specification had the effect of introducing a false moment in 
elements adjacent to the polar cap. This false moment is visible in 
Figure 6. It is important to recognize, however, that the detrimental 
effect of the erroneous specification is localized in the near-pole 
elements. 

Results indicate that the bending predictions are quite accurate 


even though a coarse mesh was used. 


C. CONCLUSIONS 

A finite element for a shell of revolution which has arbitrarily curved 
meridians has been developed and tested. Element description was based 
on displacement functions which met all requirements to guarantee mono- 
tonic convergence to an exact solution by mesh size reduction. 

Results for classical problems obtained by using the new element 
indicated that good approximations for bending stresses could be obtained 
with a very coarse mesh. Mesh size reduction will be necessary to judge 
element performance in problems where membrane action is predominant. 

The geometrical deficiency where (¢ = 0) is a limitation which must 
be recognized by users. Its effect is problem dependent. For the case 
of the spherical cap it was seen that the detrimental effect of artificially 


imposed boundary conditions at this singular point did not propagate 
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throughout the solution. On the other hand, the author was unable to 
successfully approximate stresses ina presSurized toroidal shell due to 
error propagated from similar artificial boundary conditions. Time 
precluded the determination of a method to efficiently handle this 
singularity. . 

Not withstanding this geometric limitation, results indicate that the 
new element (1) is accurate, (2) requires low computational effort and 
(3) does not require an analytical description of shell geometry. Based 
on these facts it is concluded that the curvilinear finite element can be 


efficiently used in many engineering applications. 


D, RECOMMENDATIONS 
It is recommended that: 
(1) the stress analysis program be expanded in order to handle 
realistic engineering problems; 
(2) a consistent gravity load vector be included in the stress 
analysis program; and 
(3) a study be conducted to determine a method of efficiently 


handling the singular point which exists where (¢ = 0) in the new element. 
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APPENDIX A NUMERICAL EXAMPLE OF THE STIFFNESS MATRIX 

In this appendix, a typical example of the stiffness matrix CSFE 
is presented. Eigenvalues, before and after rigid body motion was 
included, and matrix K22 are also presented. Double precision figures 


were truncated to conserve space, 
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APPENDIX B COMPUTER OUTPUT 
In this appendix, a complete set of results for the spherical cap 
analyzed in this study are presented. Results for the spherical cap were 


obtained using a 4x4 mesh. The mesh layout accompanies the problem. 
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APPENDIX C COMPUTER PROGRAM 
In this appendix the stress-strain analysis program is presented. 
User and functional information are presented through the use of comment 
cards located throughout the program. 
The program input data deck is illustrated below. Data values 


presented were those used for the spherical cap (2x2 mesh). 
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